<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><script type="text/javascript" src="https://cdn.jsdelivr.net/gh/opencobra/cobratoolbox@ffa0229fc0c01c9236bb7e961f65712443277719/latest/_static/js/iframeResizer.contentWindow.min.js"></script><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>Relaxed Flux Balance Analysis: Toy model</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: normal; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S2 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S3 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: center;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S4 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S5 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S6 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S7 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px;  }
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsErrorElement {    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsWarningElement{    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
/* Copyright 2015-2019 The MathWorks, Inc. *//* In this file, styles are not scoped to rtcContainer since they could be in the Dojo Tooltip */.diagnosticMessage-wrapper {    font-family: Menlo, Monaco, Consolas, "Courier New", monospace;    font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {    color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {    color: rgb(255,100,0);    text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {    color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {    color: rgb(230,0,0);    text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart,.diagnosticMessage-wrapper .diagnosticMessage-causePart {    white-space: normal;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {    white-space: normal;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {    white-space: normal;    word-wrap:  initial;    min-height: 18px;    max-height: 250px;    overflow: auto;}
.textElement,.rtcDataTipElement .textElement {    padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {    min-height: 16px;}
.rightPaneElement .textElement {    padding-top: 2px;    padding-left: 9px;}
.S8 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S10 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S11 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S12 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S13 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S14 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span>Relaxed Flux Balance Analysis: Toy model</span></h1><h2  class = 'S1'><span>Author: Ronan Fleming, Systems Biochemistry Group, University of Luxembourg.</span></h2><h2  class = 'S1'><span>Reviewer:</span></h2><h2  class = 'S1'><span>Introduction</span></h2><div  class = 'S2'><span>We consider a biochemical network of </span><span> </span><span>m</span><span> </span><span> molecular species and </span><span> </span><span>n</span><span> </span><span> biochemical reactions. The biochemical network is mathematically represented by a stoichiometric matrix </span><span texencoding="S\in\mathcal{Z}^{m\times n}" style="vertical-align:-5px"><img src="" width="63.5" height="19" /></span><span>. In standard notation, flux balance analysis (FBA) is the linear optimisation problem</span></div><div  class = 'S3'><span texencoding="\begin{array}{ll}
\min\limits _{v} &amp; \rho(v)\equiv c^{T}v\\
\text{s.t.} &amp; Sv=b,\\
 &amp; l\leq v\leq u,
\end{array}" style="vertical-align:-30px"><img src="" width="103.5" height="72" /></span></div><div  class = 'S2'><span>where </span><span texencoding="$c\in\Re^{n}$" style="vertical-align:-5px"><img src="" width="46.5" height="19" /></span><span> is a parameter vector that linearly combines one or more reaction fluxes to form what is termed the objective function,  and where a </span><span texencoding="$b_{i}&lt;0$" style="vertical-align:-6px"><img src="" width="39.5" height="20" /></span><span>, or  </span><span texencoding="$b_{i}&gt;0$" style="vertical-align:-6px"><img src="" width="39.5" height="20" /></span><span>, represents some fixed output, or input, of the ith molecular species. </span></div><div  class = 'S2'><span>Every FBA solution must satisfy the constraints, independent of any objective chosen to optimise over the set of constraints. It may occur that the constraints on the FBA problem are not all simultaneously feasible, i.e., the system of inequalities is infeasible. This situation might be caused by an incorrectly specified reaction bound or the absence of a reaction from the stoichiometric matrix, such that a nonzero </span><span texencoding="b\notin\mathcal{R}(S)" style="vertical-align:-5px"><img src="" width="63.5" height="19" /></span><span>. To resolve the infeasiblility, we consider a cardinality optimisation problem that seeks to minimise the number of bounds to relax, the number of fixed outputs to relax, the number of fixed inputs to relax, or a combination of all three, in order to render the problem feasible. The cardinality optimisation problem, termed </span><span style=' font-style: italic;'>relaxed flux balance analysis, </span><span>is</span></div><div  class = 'S3'><span texencoding="\begin{array}{ll}
\min\limits _{v,r,p,q} &amp; \lambda\Vert r\Vert_{0}+\alpha\Vert p\Vert_{0}+\alpha\Vert q\Vert_{0}\\
\text{s.t.} &amp; Sv+r=b\\
 &amp; l-p\leq v\leq u+q\\
 &amp; p,q,r\geq0
\end{array}" style="vertical-align:-40px"><img src="" width="183.5" height="92" /></span></div><div  class = 'S2'><span></span></div><div  class = 'S2'><span> where </span><span texencoding="$p,q\in\mathcal{R}^{n}$" style="vertical-align:-5px"><img src="" width="63.5" height="19" /></span><span> denote the relaxations of the lower and upper bounds on reaction rates of the reaction rates vector </span><span> </span><span>v, and where </span><span texencoding="$r\in\mathcal{R}^{m}$" style="vertical-align:-5px"><img src="" width="49.5" height="19" /></span><span> denotes a relaxation of the mass balance constraint. Non-negative scalar parameters </span><span> </span><span> λ </span><span> </span><span> and </span><span> </span><span> </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mi&gt;&amp;alpha;&lt;/mi&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="14.5" height="18" /></span><span> can be used to trade off between relaxation of mass balance or bound constraints. A non-negative vector parameter </span><span> </span><span> λ </span><span> </span><span> can be used to prioritise relaxation of one mass balance constraint over another, e.g, to avoid relaxation of a mass balance constraint on a metabolite that is not desired to be exchanged across the boundary of the system. A non-negative vector parameter </span><span> </span><span> </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mi&gt;&amp;alpha;&lt;/mi&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="14.5" height="18" /></span><span> </span><span> may be used to prioritise relaxation of bounds on some reactions rather than others, e.g., relaxation of bounds on exchange reactions rather than internal reactions. The optimal choice of parameters depends heavily on the biochemical context. A relaxation of the minimum number of constraints is desirable because ideally one should be able to justify the choice of bounds or choice of metabolites to be exchanged across the boundary of the system by recourse to experimental literature. This task is magnified by the number of constraints proposed to be relaxed.</span></div><h2  class = 'S1'><span>PROCEDURE: RelaxedFBA applied to a toy model</span></h2><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >clear;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >rxnForms = {</span><span style="color: rgb(170, 4, 249);">' -&gt; A'</span><span >,</span><span style="color: rgb(170, 4, 249);">'A -&gt; B'</span><span >,</span><span style="color: rgb(170, 4, 249);">'B -&gt; C'</span><span >, </span><span style="color: rgb(170, 4, 249);">'B -&gt; D'</span><span >,</span><span style="color: rgb(170, 4, 249);">'D -&gt; C'</span><span >,</span><span style="color: rgb(170, 4, 249);">'C -&gt;'</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >rxnNames = {</span><span style="color: rgb(170, 4, 249);">'R1'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R2'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R3'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R4'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R5'</span><span >, </span><span style="color: rgb(170, 4, 249);">'R6'</span><span >};</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: normal"><span >model = createModel(rxnNames, rxnNames,rxnForms);</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="D46AA51B" data-testid="output_0" data-width="420" data-height="143" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Metabolite A[c] not in model - added to the model
R1		&lt;=&gt;	A[c] 
Metabolite B[c] not in model - added to the model
R2	A[c] 	&lt;=&gt;	B[c] 
Metabolite C[c] not in model - added to the model
R3	B[c] 	&lt;=&gt;	C[c] 
Metabolite D[c] not in model - added to the model
R4	B[c] 	&lt;=&gt;	D[c] 
R5	D[c] 	&lt;=&gt;	C[c] 
R6	C[c] 	&lt;=&gt;	</div></div></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >model.lb(3) = 2;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >model.lb(4) = 2;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >model.ub(6) = 3;</span></span></div></div></div><div  class = 'S2'><span>Print the constraints</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S10'><span style="white-space: normal"><span >printConstraints(model, -1001, 1001)</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="762EB694" data-testid="output_1" data-width="420" data-height="199" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">MinConstraints:
R1	-1000
R2	-1000
R3	2
R4	2
R5	-1000
R6	-1000
maxConstraints:
R1	1000
R2	1000
R3	1000
R4	1000
R5	1000
R6	3</div></div></div></div></div><div  class = 'S2'><span>Identify the exchange reactions and biomass reaction(s) heuristically</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >model = findSExRxnInd(model,size(model.S,1),0);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%    relaxOption:    Structure containing the relaxation options:</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * internalRelax:</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 0 = do not allow to relax bounds on internal reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 1 = do not allow to relax bounds on internal reactions with finite bounds</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 2 = allow to relax bounds on all internal reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >relaxOption.internalRelax = 2;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * exchangeRelax:</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 0 = do not allow to relax bounds on exchange reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 1 = do not allow to relax bounds on exchange reactions of the type [0,0]</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 2 = allow to relax bounds on all exchange reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >relaxOption.exchangeRelax = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * steadyStateRelax:</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 0 = do not allow to relax the steady state constraint S*v = b</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 1 = allow to relax the steady state constraint S*v = b</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >relaxOption.steadyStateRelax = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * toBeUnblockedReactions - n x 1 vector indicating the reactions to be unblocked (optional)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * toBeUnblockedReactions(i) = 1 : impose v(i) to be positive</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * toBeUnblockedReactions(i) = -1 : impose v(i) to be negative</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * toBeUnblockedReactions(i) = 0 : do not add any constraint</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * excludedReactions - n x 1 bool vector indicating the reactions to be excluded from relaxation (optional)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * excludedReactions(i) = false : allow to relax bounds on reaction i</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * excludedReactions(i) = true : do not allow to relax bounds on reaction i</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * excludedMetabolites - m x 1 bool vector indicating the metabolites to be excluded from relaxation (optional)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * excludedMetabolites(i) = false : allow to relax steady state constraint on metabolite i</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * excludedMetabolites(i) = true : do not allow to relax steady state constraint on metabolite i</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div></div><div  class = 'S2'><span>Set the tolerance to distinguish between zero and non-zero flux</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(2, 128, 9);">%feasTol = getCobraSolverParams('LP', 'feasTol');</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(2, 128, 9);">%relaxOption.epsilon = feasTol/100;%*100;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    relaxOption.epsilon = 10e-6;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    relaxOption.nbMaxIteration = 1000;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    relaxOption.epsilon = 10e-6;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    relaxOption.gamma0  = 10;   </span><span style="color: rgb(2, 128, 9);">%trade-off parameter of l0 part of v</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    relaxOption.gamma1  = 1;    </span><span style="color: rgb(2, 128, 9);">%trade-off parameter of l1 part of v</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    relaxOption.lambda0 = 10;   </span><span style="color: rgb(2, 128, 9);">%trade-off parameter of l0 part of r</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    relaxOption.lambda1 = 0;    </span><span style="color: rgb(2, 128, 9);">%trade-off parameter of l1 part of r</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    relaxOption.alpha0  = 0;    </span><span style="color: rgb(2, 128, 9);">%trade-off parameter of l0 part of p and q</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    relaxOption.alpha1  = 0;     </span><span style="color: rgb(2, 128, 9);">%trade-off parameter of l1 part of p and q</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    relaxOption.theta   = 2;    </span><span style="color: rgb(2, 128, 9);">%parameter of capped l1 approximation</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div  class = 'S2'><span>Check if the model is feasible</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >FBAsolution = optimizeCbModel(model,</span><span style="color: rgb(170, 4, 249);">'max'</span><span >, 0, true);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >FBAsolution.stat == 1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    disp(</span><span style="color: rgb(170, 4, 249);">'Model is feasible. Nothing to do.'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">return</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    disp(</span><span style="color: rgb(170, 4, 249);">'Model is infeasible'</span><span >);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="069027E9" data-testid="output_2" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Model is infeasible</div></div></div></div></div><div  class = 'S11'><span>Call the relaxedFBA function, deal the solution, and set small values to zero</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S10'><span style="white-space: normal"><span >relaxOption</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="C55BB4C1" data-testid="output_3" data-width="420" data-height="76" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">relaxOption = </span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">       internalRelax: 2
       exchangeRelax: 0
    steadyStateRelax: 0
             epsilon: 1.0000e-05</div></div></div></div></div><div class="inlineWrapper"><div  class = 'S8'></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >tic;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >solution = relaxedFBA(model,relaxOption);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >timeTaken=toc;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >[v,r,p,q] = deal(solution.v,solution.r,solution.p,solution.q);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    p(p&lt;relaxOption.epsilon) = 0;</span><span style="color: rgb(2, 128, 9);">%lower bound relaxation</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    q(q&lt;relaxOption.epsilon) = 0;</span><span style="color: rgb(2, 128, 9);">%upper bound relaxation</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    r(r&lt;relaxOption.epsilon) = 0;</span><span style="color: rgb(2, 128, 9);">%steady state constraint relaxation</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div  class = 'S2'><span>The output is a solution structure with a 'stat' field reporting the solver status and a set of fields matching the relaxation of constraints given in the mathematical formulation of the relaxed flux balance problem above.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">% OUTPUT:</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%    solution:       Structure containing the following fields:</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * stat - status</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 1  = Solution found</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 0  = Infeasible</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * -1 = Invalid input</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * r - relaxation on steady state constraints S*v = b</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * p - relaxation on lower bound of reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * q - relaxation on upper bound of reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * v - reaction rate</span></span></div></div></div><div  class = 'S2'><span>Display the proposed relaxation solution</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S10'><span style="white-space: normal"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'Relaxation of steady state constraints:'</span><span >)</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="B74C29ED" data-testid="output_4" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Relaxation of steady state constraints:</div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: normal"><span >disp(r)</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="07E75E28" data-testid="output_5" data-width="420" data-height="59" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">     0
     0
     0
     0</div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: normal"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'Relaxation on lower bound of reactions:'</span><span >)</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="728509DC" data-testid="output_6" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Relaxation on lower bound of reactions:</div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: normal"><span >disp(p)</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="0DF50D40" data-testid="output_7" data-width="420" data-height="87" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">           0
           0
        1002
        1002
           0
           0</div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: normal"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,</span><span style="color: rgb(170, 4, 249);">'Relaxation on upper bound of reactions:'</span><span >)</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="E7505262" data-testid="output_8" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Relaxation on upper bound of reactions:</div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: normal"><span >disp(q)</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="C109BB0B" data-testid="output_9" data-width="420" data-height="87" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">     0
     0
     0
     0
     0
     0</div></div></div></div></div><div  class = 'S11'><span>Summarise the proposed relaxation solution</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >solution.stat == 1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >   </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >   dispCutoff=relaxOption.epsilon;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,[</span><span style="color: rgb(170, 4, 249);">'Relaxed flux balance analysis problem solved in ' </span><span >num2str(timeTaken) </span><span style="color: rgb(170, 4, 249);">' seconds.'</span><span >])</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(r),</span><span style="color: rgb(170, 4, 249);">' steady state constraints relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >     </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(p)&gt;dispCutoff &amp; ~abs(q)&gt;dispCutoff &amp; model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' internal lower bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(q)&gt;dispCutoff &amp; ~abs(p)&gt;dispCutoff &amp; model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' internal upper bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(p)&gt;dispCutoff &amp; abs(q)&gt;dispCutoff &amp; model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' internal lower and upper bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(p)&gt;dispCutoff &amp; ~abs(q)&gt;dispCutoff &amp; ~model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' external lower bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(q)&gt;dispCutoff &amp; ~abs(p)&gt;dispCutoff &amp; ~model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' external upper bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(p)&gt;dispCutoff &amp; abs(q)&gt;dispCutoff &amp; ~model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' external lower and upper bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(p)&gt;dispCutoff | abs(q)&gt;dispCutoff &amp; ~model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' external lower or upper bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    maxUB = max(max(model.ub),-min(model.lb));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    minLB = min(-max(model.ub),min(model.lb));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    intRxnFiniteBound = ((model.ub &lt; maxUB) &amp; (model.lb &gt; minLB));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(p)&gt;dispCutoff &amp; intRxnFiniteBound),</span><span style="color: rgb(170, 4, 249);">' finite lower bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(q)&gt;dispCutoff &amp; intRxnFiniteBound),</span><span style="color: rgb(170, 4, 249);">' finite upper bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    exRxn00 = ((model.ub == 0) &amp; (model.lb == 0));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(p)&gt;dispCutoff &amp; exRxn00),</span><span style="color: rgb(170, 4, 249);">' lower bounds relaxed on fixed reactions (lb=ub=0)'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(q)&gt;dispCutoff &amp; exRxn00),</span><span style="color: rgb(170, 4, 249);">' upper bounds relaxed on fixed reactions (lb=ub=0)'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    disp(</span><span style="color: rgb(170, 4, 249);">'relaxedFBA problem infeasible, check relaxOption fields'</span><span >);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="9C0A4982" data-testid="output_10" data-width="420" data-height="18" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Relaxed flux balance analysis problem solved in 0.048242 seconds.</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="D40E7BB6" data-testid="output_11" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">0 steady state constraints relaxed</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="68EFC17C" data-testid="output_12" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">2 internal lower bounds relaxed</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="3E2AFF0B" data-testid="output_13" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">0 internal upper bounds relaxed</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="49160F9F" data-testid="output_14" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">0 internal lower and upper bounds relaxed</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="519D12CE" data-testid="output_15" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">0 external lower bounds relaxed</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="A3DD2A8C" data-testid="output_16" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">0 external upper bounds relaxed</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="B0B57DC1" data-testid="output_17" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">0 external lower and upper bounds relaxed</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="32670488" data-testid="output_18" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">2 external lower or upper bounds relaxed</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="38F593A3" data-testid="output_19" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">0 finite lower bounds relaxed</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="6E9B8336" data-testid="output_20" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">0 finite upper bounds relaxed</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="A424D42F" data-testid="output_21" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">0 lower bounds relaxed on fixed reactions (lb=ub=0)</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="AD0B8427" data-testid="output_22" data-width="420" data-height="18" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">0 upper bounds relaxed on fixed reactions (lb=ub=0)</div></div></div></div><div class="inlineWrapper"><div  class = 'S13'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">return</span></span></div></div></div><h2  class = 'S1'><span>Another example</span></h2><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >rxnForms = {</span><span style="color: rgb(170, 4, 249);">' -&gt; A'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(170, 4, 249);">'A -&gt; B'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(170, 4, 249);">'A -&gt; C'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(170, 4, 249);">'D -&gt; B'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(170, 4, 249);">'E -&gt; C'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(170, 4, 249);">'E -&gt; D'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(170, 4, 249);">'E -&gt; '</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(170, 4, 249);">'A -&gt; F'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(170, 4, 249);">'G -&gt; F'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(170, 4, 249);">'H -&gt; G'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >            </span><span style="color: rgb(170, 4, 249);">'E -&gt; H'</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >rxnNames = {</span><span style="color: rgb(170, 4, 249);">'R1'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R2'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R3'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R4'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R5'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R6'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R7'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R8'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R9'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R10'</span><span >,</span><span style="color: rgb(170, 4, 249);">'R11'</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >model = createModel(rxnNames, rxnNames,rxnForms);</span></span></div></div></div><div  class = 'S2'><span>Assume all reactions are irreversible</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >model.lb(:) = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >model.ub(:) = 10;</span></span></div></div></div><div  class = 'S2'><span>Reaction R7 with bounds 1 &lt;= v_7 &lt;= 10</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: normal"><span >model.lb(7) = 1;</span></span></div></div></div><div  class = 'S2'><span>Print the constraints</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S14'><span style="white-space: normal"><span >printConstraints(model, -1001, 1001)</span></span></div></div></div><div  class = 'S2'><span>Identify the exchange reactions and biomass reaction(s) heuristically</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >model = findSExRxnInd(model,size(model.S,1),0);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%    relaxOption:    Structure containing the relaxation options:</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * internalRelax:</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 0 = do not allow to relax bounds on internal reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 1 = do not allow to relax bounds on internal reactions with finite bounds</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 2 = allow to relax bounds on all internal reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >relaxOption.internalRelax = 2;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * exchangeRelax:</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 0 = do not allow to relax bounds on exchange reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 1 = do not allow to relax bounds on exchange reactions of the type [0,0]</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 2 = allow to relax bounds on all exchange reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >relaxOption.exchangeRelax = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * steadyStateRelax:</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 0 = do not allow to relax the steady state constraint S*v = b</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 1 = allow to relax the steady state constraint S*v = b</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >relaxOption.steadyStateRelax = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * toBeUnblockedReactions - n x 1 vector indicating the reactions to be unblocked (optional)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * toBeUnblockedReactions(i) = 1 : impose v(i) to be positive</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * toBeUnblockedReactions(i) = -1 : impose v(i) to be negative</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * toBeUnblockedReactions(i) = 0 : do not add any constraint</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * excludedReactions - n x 1 bool vector indicating the reactions to be excluded from relaxation (optional)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * excludedReactions(i) = false : allow to relax bounds on reaction i</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * excludedReactions(i) = true : do not allow to relax bounds on reaction i</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * excludedMetabolites - m x 1 bool vector indicating the metabolites to be excluded from relaxation (optional)</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * excludedMetabolites(i) = false : allow to relax steady state constraint on metabolite i</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * excludedMetabolites(i) = true : do not allow to relax steady state constraint on metabolite i</span></span></div></div></div><div  class = 'S2'><span>Set the tolerance to distinguish between zero and non-zero flux</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >feasTol = getCobraSolverParams(</span><span style="color: rgb(170, 4, 249);">'LP'</span><span >, </span><span style="color: rgb(170, 4, 249);">'feasTol'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >relaxOption.epsilon = feasTol/100;</span><span style="color: rgb(2, 128, 9);">%*100;</span></span></div></div></div><div  class = 'S2'><span>Call the relaxedFBA function, deal the solution, and set small values to zero</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >tic;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >solution = relaxedFBA(model,relaxOption);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >timeTaken=toc;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >[v,r,p,q] = deal(solution.v,solution.r,solution.p,solution.q);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >0</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    p(p&lt;relaxOption.epsilon) = 0;</span><span style="color: rgb(2, 128, 9);">%lower bound relaxation</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    q(q&lt;relaxOption.epsilon) = 0;</span><span style="color: rgb(2, 128, 9);">%upper bound relaxation</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    r(r&lt;relaxOption.epsilon) = 0;</span><span style="color: rgb(2, 128, 9);">%steady state constraint relaxation</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div  class = 'S2'><span>The output is a solution structure with a 'stat' field reporting the solver status and a set of fields matching the relaxation of constraints given in the mathematical formulation of the relaxed flux balance problem above.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">% OUTPUT:</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%    solution:       Structure containing the following fields:</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * stat - status</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 1  = Solution found</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * 0  = Infeasible</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                        * -1 = Invalid input</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * r - relaxation on steady state constraints S*v = b</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * p - relaxation on lower bound of reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * q - relaxation on upper bound of reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%                      * v - reaction rate</span></span></div></div></div><div  class = 'S2'><span>Summarise the proposed relaxation solution</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >solution.stat == 1</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >   </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >   dispCutoff=relaxOption.epsilon;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >,[</span><span style="color: rgb(170, 4, 249);">'Relaxed flux balance analysis problem solved in ' </span><span >num2str(timeTaken) </span><span style="color: rgb(170, 4, 249);">' seconds.'</span><span >])</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(r),</span><span style="color: rgb(170, 4, 249);">' steady state constraints relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >     </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(p)&gt;dispCutoff &amp; ~abs(q)&gt;dispCutoff &amp; model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' internal only lower bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(q)&gt;dispCutoff &amp; ~abs(p)&gt;dispCutoff &amp; model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' internal only upper bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(p)&gt;dispCutoff &amp; abs(q)&gt;dispCutoff &amp; model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' internal lower and upper bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(p)&gt;dispCutoff &amp; ~abs(q)&gt;dispCutoff &amp; ~model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' external only lower bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(q)&gt;dispCutoff &amp; ~abs(p)&gt;dispCutoff &amp; ~model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' external only upper bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(p)&gt;dispCutoff &amp; abs(q)&gt;dispCutoff &amp; ~model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' external lower and upper bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(p)&gt;dispCutoff | abs(q)&gt;dispCutoff &amp; ~model.SIntRxnBool),</span><span style="color: rgb(170, 4, 249);">' external lower or upper bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    maxUB = max(max(model.ub),-min(model.lb));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    minLB = min(-max(model.ub),min(model.lb));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    intRxnFiniteBound = ((model.ub &lt; maxUB) &amp; (model.lb &gt; minLB));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(p)&gt;dispCutoff &amp; intRxnFiniteBound),</span><span style="color: rgb(170, 4, 249);">' finite lower bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(q)&gt;dispCutoff &amp; intRxnFiniteBound),</span><span style="color: rgb(170, 4, 249);">' finite upper bounds relaxed'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    exRxn00 = ((model.ub == 0) &amp; (model.lb == 0));</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(p)&gt;dispCutoff &amp; exRxn00),</span><span style="color: rgb(170, 4, 249);">' lower bounds relaxed on fixed reactions (lb=ub=0)'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    fprintf(</span><span style="color: rgb(170, 4, 249);">'%u%s\n'</span><span >,nnz(abs(q)&gt;dispCutoff &amp; exRxn00),</span><span style="color: rgb(170, 4, 249);">' upper bounds relaxed on fixed reactions (lb=ub=0)'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    </span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >    disp(</span><span style="color: rgb(170, 4, 249);">'relaxedFBA problem infeasible, check relaxOption fields'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><h2  class = 'S1'><span>REFERENCES</span></h2><div  class = 'S2'><span>Fleming, R.M.T., et al., Cardinality optimisation in constraint-based modelling: Application to Recon 3D (submitted), 2017.</span></div><div  class = 'S2'><span>Brunk, E. et al. Recon 3D: A resource enabling a three-dimensional view of gene variation in human metabolism. (submitted) 2017.</span></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% Relaxed Flux Balance Analysis: Toy model
%% Author: Ronan Fleming, Systems Biochemistry Group, University of Luxembourg.
%% Reviewer:
%% Introduction
% We consider a biochemical network of  m  molecular species and  n  biochemical 
% reactions. The biochemical network is mathematically represented by a stoichiometric 
% matrix $S\in\mathcal{Z}^{m\times n}$. In standard notation, flux balance analysis 
% (FBA) is the linear optimisation problem
% 
% $$\begin{array}{ll}\min\limits _{v} & \rho(v)\equiv c^{T}v\\\text{s.t.} & 
% Sv=b,\\ & l\leq v\leq u,\end{array}$$
% 
% where $$c\in\Re^{n}$$ is a parameter vector that linearly combines one or 
% more reaction fluxes to form what is termed the objective function,  and where 
% a $$b_{i}<0$$, or  $$b_{i}>0$$, represents some fixed output, or input, of the 
% ith molecular species. 
% 
% Every FBA solution must satisfy the constraints, independent of any objective 
% chosen to optimise over the set of constraints. It may occur that the constraints 
% on the FBA problem are not all simultaneously feasible, i.e., the system of 
% inequalities is infeasible. This situation might be caused by an incorrectly 
% specified reaction bound or the absence of a reaction from the stoichiometric 
% matrix, such that a nonzero $b\notin\mathcal{R}(S)$. To resolve the infeasiblility, 
% we consider a cardinality optimisation problem that seeks to minimise the number 
% of bounds to relax, the number of fixed outputs to relax, the number of fixed 
% inputs to relax, or a combination of all three, in order to render the problem 
% feasible. The cardinality optimisation problem, termed _relaxed flux balance 
% analysis,_ is
% 
% $$\begin{array}{ll}\min\limits _{v,r,p,q} & \lambda\Vert r\Vert_{0}+\alpha\Vert 
% p\Vert_{0}+\alpha\Vert q\Vert_{0}\\\text{s.t.} & Sv+r=b\\ & l-p\leq v\leq u+q\\ 
% & p,q,r\geq0\end{array}$$
% 
% 
% 
% where $$p,q\in\mathcal{R}^{n}$$ denote the relaxations of the lower and upper 
% bounds on reaction rates of the reaction rates vector  v, and where $$r\in\mathcal{R}^{m}$$ 
% denotes a relaxation of the mass balance constraint. Non-negative scalar parameters   
% λ   and   $\alpha \;$ can be used to trade off between relaxation of mass balance 
% or bound constraints. A non-negative vector parameter   λ   can be used to prioritise 
% relaxation of one mass balance constraint over another, e.g, to avoid relaxation 
% of a mass balance constraint on a metabolite that is not desired to be exchanged 
% across the boundary of the system. A non-negative vector parameter   $\alpha 
% \;$  may be used to prioritise relaxation of bounds on some reactions rather 
% than others, e.g., relaxation of bounds on exchange reactions rather than internal 
% reactions. The optimal choice of parameters depends heavily on the biochemical 
% context. A relaxation of the minimum number of constraints is desirable because 
% ideally one should be able to justify the choice of bounds or choice of metabolites 
% to be exchanged across the boundary of the system by recourse to experimental 
% literature. This task is magnified by the number of constraints proposed to 
% be relaxed.
%% PROCEDURE: RelaxedFBA applied to a toy model

clear;
rxnForms = {' -> A','A -> B','B -> C', 'B -> D','D -> C','C ->'};
rxnNames = {'R1','R2','R3','R4','R5', 'R6'};
model = createModel(rxnNames, rxnNames,rxnForms);
model.lb(3) = 2;
model.lb(4) = 2;
model.ub(6) = 3;
%% 
% Print the constraints

printConstraints(model, -1001, 1001)
%% 
% Identify the exchange reactions and biomass reaction(s) heuristically

model = findSExRxnInd(model,size(model.S,1),0);
%    relaxOption:    Structure containing the relaxation options:
%
%                      * internalRelax:
%
%                        * 0 = do not allow to relax bounds on internal reactions
%                        * 1 = do not allow to relax bounds on internal reactions with finite bounds
%                        * 2 = allow to relax bounds on all internal reactions
relaxOption.internalRelax = 2;
%
%                      * exchangeRelax:
%
%                        * 0 = do not allow to relax bounds on exchange reactions
%                        * 1 = do not allow to relax bounds on exchange reactions of the type [0,0]
%                        * 2 = allow to relax bounds on all exchange reactions
relaxOption.exchangeRelax = 0;
%                      * steadyStateRelax:
%
%                        * 0 = do not allow to relax the steady state constraint S*v = b
%                        * 1 = allow to relax the steady state constraint S*v = b
relaxOption.steadyStateRelax = 0;
%                      * toBeUnblockedReactions - n x 1 vector indicating the reactions to be unblocked (optional)
%
%                        * toBeUnblockedReactions(i) = 1 : impose v(i) to be positive
%                        * toBeUnblockedReactions(i) = -1 : impose v(i) to be negative
%                        * toBeUnblockedReactions(i) = 0 : do not add any constraint
%
%                      * excludedReactions - n x 1 bool vector indicating the reactions to be excluded from relaxation (optional)
%
%                        * excludedReactions(i) = false : allow to relax bounds on reaction i
%                        * excludedReactions(i) = true : do not allow to relax bounds on reaction i
%
%                      * excludedMetabolites - m x 1 bool vector indicating the metabolites to be excluded from relaxation (optional)
%
%                        * excludedMetabolites(i) = false : allow to relax steady state constraint on metabolite i
%                        * excludedMetabolites(i) = true : do not allow to relax steady state constraint on metabolite i
%
%% 
% Set the tolerance to distinguish between zero and non-zero flux

if 1
    %feasTol = getCobraSolverParams('LP', 'feasTol');
    %relaxOption.epsilon = feasTol/100;%*100;
    relaxOption.epsilon = 10e-6;
else
    relaxOption.nbMaxIteration = 1000;
    relaxOption.epsilon = 10e-6;
    relaxOption.gamma0  = 10;   %trade-off parameter of l0 part of v
    relaxOption.gamma1  = 1;    %trade-off parameter of l1 part of v
    relaxOption.lambda0 = 10;   %trade-off parameter of l0 part of r
    relaxOption.lambda1 = 0;    %trade-off parameter of l1 part of r
    relaxOption.alpha0  = 0;    %trade-off parameter of l0 part of p and q
    relaxOption.alpha1  = 0;     %trade-off parameter of l1 part of p and q
    relaxOption.theta   = 2;    %parameter of capped l1 approximation
end
%% 
% Check if the model is feasible

FBAsolution = optimizeCbModel(model,'max', 0, true);
if FBAsolution.stat == 1
    disp('Model is feasible. Nothing to do.');
    return
else
    disp('Model is infeasible');
end
%% 
% Call the relaxedFBA function, deal the solution, and set small values to zero

relaxOption

tic;
solution = relaxedFBA(model,relaxOption);
timeTaken=toc;
[v,r,p,q] = deal(solution.v,solution.r,solution.p,solution.q);
if 0
    p(p<relaxOption.epsilon) = 0;%lower bound relaxation
    q(q<relaxOption.epsilon) = 0;%upper bound relaxation
    r(r<relaxOption.epsilon) = 0;%steady state constraint relaxation
end
%% 
% The output is a solution structure with a 'stat' field reporting the solver 
% status and a set of fields matching the relaxation of constraints given in the 
% mathematical formulation of the relaxed flux balance problem above.

% OUTPUT:
%    solution:       Structure containing the following fields:
%                      * stat - status
%                        * 1  = Solution found
%                        * 0  = Infeasible
%                        * -1 = Invalid input
%                      * r - relaxation on steady state constraints S*v = b
%                      * p - relaxation on lower bound of reactions
%                      * q - relaxation on upper bound of reactions
%                      * v - reaction rate
%% 
% Display the proposed relaxation solution

fprintf('%s\n','Relaxation of steady state constraints:')
disp(r)
fprintf('%s\n','Relaxation on lower bound of reactions:')
disp(p)
fprintf('%s\n','Relaxation on upper bound of reactions:')
disp(q)
%% 
% Summarise the proposed relaxation solution

if solution.stat == 1
   
   dispCutoff=relaxOption.epsilon;

    fprintf('%s\n',['Relaxed flux balance analysis problem solved in ' num2str(timeTaken) ' seconds.'])
    
    fprintf('%u%s\n',nnz(r),' steady state constraints relaxed');
     
    fprintf('%u%s\n',nnz(abs(p)>dispCutoff & ~abs(q)>dispCutoff & model.SIntRxnBool),' internal lower bounds relaxed');
    fprintf('%u%s\n',nnz(abs(q)>dispCutoff & ~abs(p)>dispCutoff & model.SIntRxnBool),' internal upper bounds relaxed');
    fprintf('%u%s\n',nnz(abs(p)>dispCutoff & abs(q)>dispCutoff & model.SIntRxnBool),' internal lower and upper bounds relaxed');

    fprintf('%u%s\n',nnz(abs(p)>dispCutoff & ~abs(q)>dispCutoff & ~model.SIntRxnBool),' external lower bounds relaxed');
    fprintf('%u%s\n',nnz(abs(q)>dispCutoff & ~abs(p)>dispCutoff & ~model.SIntRxnBool),' external upper bounds relaxed');
    fprintf('%u%s\n',nnz(abs(p)>dispCutoff & abs(q)>dispCutoff & ~model.SIntRxnBool),' external lower and upper bounds relaxed');
    
    fprintf('%u%s\n',nnz(abs(p)>dispCutoff | abs(q)>dispCutoff & ~model.SIntRxnBool),' external lower or upper bounds relaxed');
    
    maxUB = max(max(model.ub),-min(model.lb));
    minLB = min(-max(model.ub),min(model.lb));
    intRxnFiniteBound = ((model.ub < maxUB) & (model.lb > minLB));
    fprintf('%u%s\n',nnz(abs(p)>dispCutoff & intRxnFiniteBound),' finite lower bounds relaxed');
    fprintf('%u%s\n',nnz(abs(q)>dispCutoff & intRxnFiniteBound),' finite upper bounds relaxed');
    
    exRxn00 = ((model.ub == 0) & (model.lb == 0));
    fprintf('%u%s\n',nnz(abs(p)>dispCutoff & exRxn00),' lower bounds relaxed on fixed reactions (lb=ub=0)');
    fprintf('%u%s\n',nnz(abs(q)>dispCutoff & exRxn00),' upper bounds relaxed on fixed reactions (lb=ub=0)');
    
else
    disp('relaxedFBA problem infeasible, check relaxOption fields');
end
return
%% Another example

rxnForms = {' -> A',...
            'A -> B',...
            'A -> C',...
            'D -> B',...
            'E -> C',...
            'E -> D',...
            'E -> ',...
            'A -> F',...
            'G -> F',...
            'H -> G',...
            'E -> H'};
rxnNames = {'R1','R2','R3','R4','R5','R6','R7','R8','R9','R10','R11'};
model = createModel(rxnNames, rxnNames,rxnForms);
%% 
% Assume all reactions are irreversible

model.lb(:) = 0;
model.ub(:) = 10;
%% 
% Reaction R7 with bounds 1 <= v_7 <= 10

model.lb(7) = 1;
%% 
% Print the constraints

printConstraints(model, -1001, 1001)
%% 
% Identify the exchange reactions and biomass reaction(s) heuristically

model = findSExRxnInd(model,size(model.S,1),0);
%    relaxOption:    Structure containing the relaxation options:
%
%                      * internalRelax:
%
%                        * 0 = do not allow to relax bounds on internal reactions
%                        * 1 = do not allow to relax bounds on internal reactions with finite bounds
%                        * 2 = allow to relax bounds on all internal reactions
relaxOption.internalRelax = 2;
%
%                      * exchangeRelax:
%
%                        * 0 = do not allow to relax bounds on exchange reactions
%                        * 1 = do not allow to relax bounds on exchange reactions of the type [0,0]
%                        * 2 = allow to relax bounds on all exchange reactions
relaxOption.exchangeRelax = 0;
%                      * steadyStateRelax:
%
%                        * 0 = do not allow to relax the steady state constraint S*v = b
%                        * 1 = allow to relax the steady state constraint S*v = b
relaxOption.steadyStateRelax = 0;
%                      * toBeUnblockedReactions - n x 1 vector indicating the reactions to be unblocked (optional)
%
%                        * toBeUnblockedReactions(i) = 1 : impose v(i) to be positive
%                        * toBeUnblockedReactions(i) = -1 : impose v(i) to be negative
%                        * toBeUnblockedReactions(i) = 0 : do not add any constraint
%
%                      * excludedReactions - n x 1 bool vector indicating the reactions to be excluded from relaxation (optional)
%
%                        * excludedReactions(i) = false : allow to relax bounds on reaction i
%                        * excludedReactions(i) = true : do not allow to relax bounds on reaction i
%
%                      * excludedMetabolites - m x 1 bool vector indicating the metabolites to be excluded from relaxation (optional)
%
%                        * excludedMetabolites(i) = false : allow to relax steady state constraint on metabolite i
%                        * excludedMetabolites(i) = true : do not allow to relax steady state constraint on metabolite i
%% 
% Set the tolerance to distinguish between zero and non-zero flux

feasTol = getCobraSolverParams('LP', 'feasTol');
relaxOption.epsilon = feasTol/100;%*100;
%% 
% Call the relaxedFBA function, deal the solution, and set small values to zero

tic;
solution = relaxedFBA(model,relaxOption);
timeTaken=toc;
[v,r,p,q] = deal(solution.v,solution.r,solution.p,solution.q);
if 0
    p(p<relaxOption.epsilon) = 0;%lower bound relaxation
    q(q<relaxOption.epsilon) = 0;%upper bound relaxation
    r(r<relaxOption.epsilon) = 0;%steady state constraint relaxation
end
%% 
% The output is a solution structure with a 'stat' field reporting the solver 
% status and a set of fields matching the relaxation of constraints given in the 
% mathematical formulation of the relaxed flux balance problem above.

% OUTPUT:
%    solution:       Structure containing the following fields:
%                      * stat - status
%                        * 1  = Solution found
%                        * 0  = Infeasible
%                        * -1 = Invalid input
%                      * r - relaxation on steady state constraints S*v = b
%                      * p - relaxation on lower bound of reactions
%                      * q - relaxation on upper bound of reactions
%                      * v - reaction rate
%% 
% Summarise the proposed relaxation solution

if solution.stat == 1
   
   dispCutoff=relaxOption.epsilon;

    fprintf('%s\n',['Relaxed flux balance analysis problem solved in ' num2str(timeTaken) ' seconds.'])
    
    fprintf('%u%s\n',nnz(r),' steady state constraints relaxed');
     
    fprintf('%u%s\n',nnz(abs(p)>dispCutoff & ~abs(q)>dispCutoff & model.SIntRxnBool),' internal only lower bounds relaxed');
    fprintf('%u%s\n',nnz(abs(q)>dispCutoff & ~abs(p)>dispCutoff & model.SIntRxnBool),' internal only upper bounds relaxed');
    fprintf('%u%s\n',nnz(abs(p)>dispCutoff & abs(q)>dispCutoff & model.SIntRxnBool),' internal lower and upper bounds relaxed');

    fprintf('%u%s\n',nnz(abs(p)>dispCutoff & ~abs(q)>dispCutoff & ~model.SIntRxnBool),' external only lower bounds relaxed');
    fprintf('%u%s\n',nnz(abs(q)>dispCutoff & ~abs(p)>dispCutoff & ~model.SIntRxnBool),' external only upper bounds relaxed');
    fprintf('%u%s\n',nnz(abs(p)>dispCutoff & abs(q)>dispCutoff & ~model.SIntRxnBool),' external lower and upper bounds relaxed');
    
    fprintf('%u%s\n',nnz(abs(p)>dispCutoff | abs(q)>dispCutoff & ~model.SIntRxnBool),' external lower or upper bounds relaxed');
    
    maxUB = max(max(model.ub),-min(model.lb));
    minLB = min(-max(model.ub),min(model.lb));
    intRxnFiniteBound = ((model.ub < maxUB) & (model.lb > minLB));
    fprintf('%u%s\n',nnz(abs(p)>dispCutoff & intRxnFiniteBound),' finite lower bounds relaxed');
    fprintf('%u%s\n',nnz(abs(q)>dispCutoff & intRxnFiniteBound),' finite upper bounds relaxed');
    
    exRxn00 = ((model.ub == 0) & (model.lb == 0));
    fprintf('%u%s\n',nnz(abs(p)>dispCutoff & exRxn00),' lower bounds relaxed on fixed reactions (lb=ub=0)');
    fprintf('%u%s\n',nnz(abs(q)>dispCutoff & exRxn00),' upper bounds relaxed on fixed reactions (lb=ub=0)');
    
else
    disp('relaxedFBA problem infeasible, check relaxOption fields');
end
%% REFERENCES
% Fleming, R.M.T., et al., Cardinality optimisation in constraint-based modelling: 
% Application to Recon 3D (submitted), 2017.
% 
% Brunk, E. et al. Recon 3D: A resource enabling a three-dimensional view of 
% gene variation in human metabolism. (submitted) 2017.
##### SOURCE END #####
-->
</div></body></html>